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ABSTRACT 


The basic ideas of the construction of numerically generated body-fitted 
coordinate systems for the finite difference solution of partial differential equa- 
tions are discussed. Algebraic grid generation and grid generation based on 
partial differential equations have been discussed in detail with specific refer- 
ence to Transfinite Interpolation, Laplace, Thomas-Middlecoff, Sorenson, and 
Modified Sorenson method. Some parameters which affect the ‘goodness’ of a 
grid are prescribed. Finally, the grid generated by various methods for various 
test cases are compared, to evaluate their performance on H-type grids gen- 
erated for viscous, internal flows. It was found that Thomas-Middlecoff gives 
the best performance with Transfinite Interpolation next best. 
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A 
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Cell area 

Arbitrary constants 
d Arbitrary constants 
Denotes distance 

Grid node corresponding to (2 — l)Aif 
Transformation Jacobian matrix 
Grid node corresponding to {j — 1)A77 
Source terms of poisson’s equation 
A vector denotes a point (x, y) 

Denotes distance 
Denotes distance 
Denotes distance 
Cartesian coordinates 

Coefficients in transformed poisson’s equation 
Angle between grid lines (deg) 

Curvilinear coordinates 
Thomas-Middlecoff grid control functions 
Grid interval in ^ direction 
Grid interval in 77 direction 


Superscripts 

Represents iteration level 
Represents computational domain 


Represents 

Represents 

Represents 

Represents 

Represents 

Represents 

Represents 

Represents 


Subscripts 
values at boundary 
interior point 
maximum value 
minimum value 

partial derivatives with respect to x keeping y constant 
partial derivatives with respect to y keeping x constant 
partial derivatives with respect to ^ keeping 77 constant 
partial derivatives with respect to 77 keeping ^ constant 



Chapter 1 


INTRODUCTION 


The generation of body-fitted coordinate systems plays an important role in the 
accurate computation of numerical solutions of partial differential equations on 
regions with arbitrary shaped boundaries. A grid which is not well suited for 
the problem can lead to unsatisfactory results. In some applications, improper 
choice of grid point locations can lead to instability or lack of convergence, 
While with a proper selection of grid points, the accuracy of the numerical 
simulation can be increased without increased computational effort. 


1.1 Classification of Grids 

Grids cam be classified into structured grids and unstructured grids. In struc- 
tured grid generation, the grid points in the comple.x physical domain have a 
one to one correspondence with the grid points of a rectangular grid of the 
same size in the computational domain. Structured grids are primarily used 
for obtaining solutions based on finite difference or finite volume methods. 
The number of nearest neighbours for each interior grid point is fixed in a 
structured grid (i.e. 2, 4, 6 for one, two and three dimensions respectively). 



In unstructured grid generation, no definite structure is maintained be- 
tween any typical node and its neighbours. Therefore, unstructured grids axe 
best suited for finite element applications. 

We are interested only in structured grid generation. So unstructured 
grid generation will not be discussed here after. 

1.2 Importance of Body-fitted Coordinate Sys- 
tems 

With a transformation approach, all computations can be done on a square 
grid in a transformed domain regardless of the shape of the physical domain. In 
rectangular enveloped grid systems, the cells adjacent to a complex boundary 
will be incomplete, as they will approximate the boundary using rectangular 
cells and interpolation is required over these cells. In body-fitted coordinate 
systems, the body surface is a boundary in the computational domain which 
enables better implementation of boundary conditions. Thus body-fitted grid 
system avoids the need for modifying the solver for a new geometry. 


1.3 Present Work 

Many packages are available for structured grid generation. Each of them 
will be based on one or two grid generation methods and are suitable only 
for certain kinds of problems. So it is difficult to get a package for universal 
application. These packages are generally a black box, and we do not exactly 
know what is going inside. Also, for a research group, it is desirable to have 
a package of our own to suit our needs which can be fine tuned as and when 
desired. 
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So, in our present work, we developed a code in FORTRAN 77 for struc- 
tured grid generation in two dimensions using various methods. The code can 
also evaluate the ‘goodness’ of the grid using various criteria which are dis- 
cussed in chapter 2. Finally, the grids generated by various techniques are 
compared in chapter 4. 



Chapter 2 


PRINCIPLES OF GRID 
GENERATION 


The fundamental problem in grid generation is to systematically choose a grid 
which allows for an accurate numerical solution of the governing equations, 
with the minimum computation. The physical domain on which the problem is 
to be solved, is often mapped onto a rectangular computational domain. Then 
the problem is solved in the computational domain by solving the governing 
equations, suitably transformed, in the rectangular grid. In multiblock grid 
generation, separate sections of the physical domain are mapped onto separate 
squares or rectangles. 

Figure 2.1 shows a structured grid in the physical and transformed do- 
main. Note that for each grid point (z, j) in the physical domain there is a 
corresponding grid point in the transformed domain. Also for each grid line 
{A — A or B — B) there ane corresponding lines (A' — A' , B' — B') in the 
transformed domain. Further, the grid lines in the transformed domain are 
lines of constant ^ or lines of constant rj. 


Various aspects of structured grid generation is discussed in this chapter. 
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(a) Physical domain Computational domain 


Figure 2.1: 


2.1 Requirements of a Computational Grid 

For better results of numerical simulations, there are some disirable features 
required of the grid point distribution. Some of them are. 


1. The mapping from physical domain to computational domain must be 
one-to-one 

2. The grid lines should be smooth so as to provide continuous transforma- 
tion parameters (which appear in the transformed governing equations) . 

3 Grid points should be closely spaced in the physical domain where large 
numerical errors are expected, (eg. fine grids are necessary in the bound- 
ary layer region). 

4. Excessive skewness should be avoided. If the angles between grid lines 
are close to zero, the numerical simulation will be highly inaccurate due 
to loss of independence between the grid lines. 

5 It is desirable to make the intersecting grid lines orthogonal in the bound- 
ary for easy implementation of boundary conditions. 
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It is not possible to meet all the demands simultaneously . A particu- 
lar property may be very important, and the requirements may change from 
problem to problem with some feasible compromise on the others The chief 
concern m body-fitted grid generation is to ensure smoothness while efficient 
grid spacing is the primary interest for adaptive grid generation. Apart from 
the above requirements, grid generation process should be automatic so as to 
avoid the drudgery of giving grid data manually. 

2.2 Classification of Structured Grids 

Two dimensional grids which map the physical domain to computational do- 
main are generally classified into three categories. They are: 

(a) 0-type grids 

(b) C-type grids 

(c) H-type grids 



* 


(a) Physical domain 


(b) Computational domain 


Figure 2.2: 
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For a given geometry any one of these configurations can be obtained by 
suitable mapping Consider the domain shown in figure 2 2(a), which contains 
a ‘hole’ (eg a solid body) which is to be excluded from the computational 
domain shown on the right 


y 
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(a) Grids in physical domain 


A B' 

(b) computational domain 


Figure 2.3: 0 - Type Grid 


In 0-type grid generation, the extenor of the physical domain is mapped 
onto one side of the computational domain. To get 0-grids for this domain, 
shown in figure 2.2(a), a branch cut is introduced in the domain and the points 
A,B,C,D are chosen in either sides of the branch cut as shown in the figure 
2.3(a). Then AB. BC, CD, DA are mapped onto A'B’, B'C, C'D\ D'A' 
respectively. Thus the two sides of the branch cut (AD, BC) are mapped 
onto two opposite sides of the computational domain, and the line A'B' and 
B'C physically correspond to the same line (branch cut) It can be seen that in 
this case, 77 -constant lines forms 0-type configuration in the physical domain. 

In C-type grid generation, the outer boundary is mapped onto three sides 
of the computational domain. To generate C-type grids for the same domain 
shown in figure 2.2(a), a branch cut is introduced and the points A, B are 
chosen on either side where the branch cut meets the outer boundary. The 
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Figure 2.4. C - Type Grid 


points C, D are suitably selected in the outer boundary and then AB, BC, CD, 
DA are mapped onto A'B' , B'C, CD', D'A' respectively. The forward sweep 
branch cut (AP), object surface (PQ) and reverse sweep branch cut (QB) form 
one side of transformed region. The solid body is reduced to the line P'Q' at 
the bottom boundary of the transformed domain. It can be seen that in this 
case, 77-constant lines form a ‘C’ figure in the physical domain. 




Figure 2.5: H - Type Grid 





2.3 Transformation Relations 
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For H-type grid generation, the outer boundary is mapped onto four 
sides of the computational domain. To generate H-grids for the same domain 
shown in figure 2.2(a), branch cuts are introduced on either sides of the object 
(figure 2.5(a)) and the upper and lower portions ABCQPD and EPQFGH are 
separately mapped onto A'B'C'Q'P'D' and E' P'Q'F'G' H' . The object in the 
physical domain reduces to a line in the computational domain. 

C-grids and O-grids are mostly used for external flow calculations such 
as flow over airfoils, cascades etc. . while H-grids are mostly used for internal 
flow calculations. 


2.3 Transformation Relations 

To construct numerical solution of partial differential equations on regions of 
arbitrary shape, the computations have to be carried out in the transformed 
region with the transformed coordinates (i^, 77) as the independent variables. 
The coordinate in the physical domain will be denoted by (x, y) and those 
of the transformed region by rj). Any point on the physical domain has a 
single corresponding point on the transformed region, and vice-versa. In this 
section, the transformation relations between cartesian coordinate system and 
body-fitted curvilinear coordinate (^, 77) system are discussed. 

The body-fitting coordinate transformation will be 

^ y) , 

V = y) > 

the bounds of the computational domain are defined by 


o<e<i, 

0 <77 < 1 . 
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From the chain rule 



Vx Vy 



(2.1) 


where and rjx, rjy are the partial derivatives of rj with respect to x, y. 

Also we have the inverse relation ship 


X = x(^, 77 ) , 

y = y(f, v) ■ 


Again using chain rule 


dx 

dy 


Xi x^ / 1 

Vi Vv \ \dn i 


substituting Equation (2.1) in Equation (2.2) 


dx 

dy 


X^ Xy! 

Vi Vn 


^x 

nx rjy 



( 2 . 2 ) 


(2.3) 


we see that the two matrices on the right must be inverses of each other, i.e., 

(2.4) 


1 



1 

Vv 

Xt) 

_ T]y _ 


1 

J 

. -y^ 

xi 


where | J |= — x^y,, is the determinant of the Jacobian matrix 


[J] = 


Xr, 

yr, 


For a one-to-one mapping, the determinant | J | should be finite and 
non-zero. It can be shown that the cell area in the physical domain 

dA =\ J \ d^drj 

while the corresponding area in the transformed region is d^drj. Thus the 
determinant | J \ represents the local area scaling factor. 
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From the Equation 2.4, The relations between the partial derivatives 
of curvihnear coordinates (if, 77 ) and that of cartesian coordinates (x, y) are 
obtained as 


Vv = 


Vv 

\J\ 

Z5il 

\J\ 

ZK 

\J\ 


and 


(2.5) 

( 2 . 6 ) 

(2.7) 


Vy = 


ur 


(2.8) 


These relations are very useful in transforming the governing equations from 
the (x, y) coordinate system to (f, rj) system. 


The governing equations will be in terms of the physical domain variables 
(x, y). However, as they will be solved in the transformed domain, these 
equations need to be put in terms of (f, rj). This is done by the mapping 


X = X(f, T]) 

y = r}) 


and 


dx , 




. dx 


A 


+ ( 1 

(d£\ ^ 
\dr)), \J\ 




etc. 


so that the equations are finally expressed only in terms of (f, 77 ) derivatives 
with respect to (f, 77 ). The transformed equations are normally much more 
complicated than the original equations. This price has to be paid for the 
simplification in geometry. 
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2.4 Grid Generation Techniques 

The techniques available for structured grid generation can be divided into 
three major categories. They are 

1. Algebraic grid generation 

2. Conformal mapping 

3. Methods using Partial Differential Equations 

2.4.1 Algebraic grid generation 

Among the grid generation techniques, algebraic grid generation computation- 
ally most efficient. These methods are based on mathematical interpolation 
functions. Explicit control of grid shape and grid spacing is also possible. The 
disadvantage of these methods is that the slope discontinuities present at the 
boundaries tend to propagate inside and hence smoothness is not assured. 

Some of the well known algebraic grid generation methods are transfinite 
interpolation, multisurface method and two boundary method. Transfinite 
interpolation is most widely used method. This method is discussed in detail 
in the next section. 

2.4.2 Conformal mapping 

Conformal mapping provides a smooth orthogonal grid except at a few singular 
points. This method can be used only for two dimensional problems. In this 
thesis, we are not concerned with this method and hence it is not discussed 
here. 
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2.4.3 Grid generation based on Partial Differential Equa- 
tions 

These methods are the most commonly used. These methods can be fur- 
ther classified into elHptic, hyperbolic and parabolic grid generation depending 
upon the generation equations used. The grids generated using this methods 
are smooth, and control of orthogonality and grid spacing is possible. But 
these methods are computationally expensive. 

Elliptic grid generation uses as grid generation equations Laplace or Pois- 
son’s equations, where all four boundaries need to be specified. Hence this 
method is used for closed geometries. Elliptic grid generation will be discussed 
in the next chapter in detail. 

In hyperbolic grid generation, the outer boundary need not be specified. 
The grid lines generated using hyperbolic equations are orthogonal near the 
boundary. These methods are used to generate grids for domains that are 
open. 


Parabolic grid generation is not frequently used. Achieving orthogonality 
is not effective as hyperbolic method. Also, control functions are not developed 
for this method. Both parabolic and hyperbolic grid generation are faster than 
elhptic grid generation. 


2.5 Transfinite Interpolation 

This is a powerful but simple algebraic technique for grid generation and is fre- 
quently used. The interior coordinates are generated by series of unidirectional 
interpolation between prescribed boundary data. 

Consider a 2-D physical domain ABCD shown in figure 2.6(a). This is 



2.5 Transfinite Interpolation 


14 


to be mapped onto a square computational domain with equally spaced grid 
points as shown in figure 2.6(e) 

First, grid points are placed on the boundary in the physical domain. 
This is done by choosing the grid points on the boundaries according to some 
sensible criterion Fig 2.6(b). The number of grid points on opposite sides 
of the domain should be equal. These boundary grid points have a one-to- 
one correspondence with the boundary points on the transformed domain Fig 
2.6(e). The (^, 77) values of each grid point in Fig 2.6(e) are known as points 
are equally spaced and the domain is a unit square in (^, rj) space. Now we 
need to find the (x, y)coordinates for every interior grid point in Fig 2.6(b), 
corresponding to each grid point in Fig 2.6(e). 

Now unidirectional interpolation is applied in any one of the direction, 
say the ^-direction. For this we do the following. For every rj = constant line 
in figure 2.6(e) interpolate the (x, y) coordinates of a point P using 

Xp = (1 — 

j/p = (1 - ^)yE + ^Vf 

Where E and F are the end of the 77 — constant line. This interpolation 
is carried out for all 77 — constant lines including the boundary lines 77 = 0 and 
77 = 1. These interpolated points are shown in the figure 2.6(c). It can be seen 
in this figure that, while the boundary points at the ^ = 0 and ^ = 1 
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(a) Physical domain 


y 




(c) After unidirectional interpolation in 
(f - direction 



X 


(d) Final grid 


F| 



Figure 2.6. Transfinite Interpolation 
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boundaries (AD and BC) the interpolated boundary points match exactly with 
the chosen boundary points, at the boundaries AB and CD, the grid points 
obtained by unidirectional interpolation need not coincide with the chosen 
boundary points and may even fall outside the domain. 

This anomaly is corrected by unidirectional interpolation in the 77 di- 
rection. These interpolation are made along each f — constant line. The 
corrections Axp and Ayp for each point P are made such that at the 77 = 0 
and 7 = 1 boundaries the corrections are exactly the difference (Axa, Aya) 
and [Axff, Ayjj) which would bring the ‘boundary’ points obtained by the 
first interpolation exactly to the chosen boundary points. 

Axp = (1 — rf)AxG + rjAxjj 

Ayp = (1 - ri)AyG + yAyn 

Where G and H are the end points of the ^ = constant line on which P lies. 
Axg, Axff, AyG and Ayp are the corrections for the boundary points G and 
H and is given by 

Axg = Xq - xg Aya = yG- 

Axff = x'fj -XH Ayp = y'fj — xp 

Where Xq - xg etc. are the values obtained by the first interpolation and xq 
etc. are the values chosen for the boundary points. 

The final values of the coordinates of P are obtained as 

Xp = x'p — Axp 

yp = y'p- ^yp 

The final mesh in the physical domain will appear as shown in the figure 2 . 6 (d), 
each point on which has a corresponding point in the transformed domain. 
Now as the (xp, yp) values in the physical domain, and the (^p, 7 p) values 
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in the transformed domain, are known for every point P, the relationship 
X — X (if, 77 ) and y = y (f , 77 ) are also known. 

The quality of the mesh generated entirely depends upon the chosen 
boundary points distribution. So care should be taken m the placement of 
the boundary points. Transfinite interpolation can be applied to 2 -D and 3-D 
geometries and for all grid configurations ie 0 - type. C - type and H type 
grids. 

2.6 Evaluation of Grid Quality 

Truncation error of numerical solutions depends upon the grid upon which the 
discretized approximation of the differential equation is solved. A desirable 
grid should give minimum truncation error. Certain properties of the grid are 
found to have a strong effects on the truncation error [5]. Parameters which 
affect these properties can be used to express the goodness of a grid. The 
parameters are 

1 . Transformation Jacobian 

2. Skewness 

3. Aspect ratio 

4. Adjacent cell ratio 

5. Nearness parameter 

Transformation Jacobian 

As we have already mentioned, the transformation Jacobian represents the 
local area scaling factor. So the Jacobian should not become zero for the 
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mapping to be one-to-one. 

Skewness 

Skewness is a measure of nonorthogonality between two intersecting grid lines. 
Higher the skewness, higher will be the truncation error due to loss of inde- 
pendence between the grid lines. Many methods are available to evaluate the 
cell skewness. They are 

(a) Intersection Angle 

The intersection angle 9 between two grid lines, as shown in figure 2.7, is 
calculated. It is desirable that the intersection angle be 90°, but this is not 
always possible. At least it is desirable to maintain the intersection angle in 
the range 

45 ° <9 < 135° 



Figure 2.7: 


(b) Based on length of diagonals 

In this method skewness is calculated by using the length of the diagonals of 
a grid cell. For a grid cell ABCD as shown in figure 2.8 skewness is calculated 
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as 


skewness = 1 — 


length of smaller diagonal BD 
length of larger diagonal AC 


D 



C 


(c) Based on length of half diagonal 


For a grid cell of shape as shown in figure 2.9, the previous method of calcu- 
lating cell skewness will fail to give true information about the cell skewness 
for the diagonals are equal and the skewness will be zero, although it can be 


seen the quadrilateral is quite skewed. For such cases, the length of the half 

diagonals is used to calculate the skewness. 

length of smallest half diagonal OD 


skewness = 1— 


length of larger half diagonal of the opposite diagonal OA 


(d) Based on area of triangles 


In this method, skewness is calculated using the area of the triangles formed 
by the diagonals and the four sides of a grid cell. 


For a grid shown in figure 2 9 

area of the smallest triangle AODC 


skewness 


area of biggest triangle AOAB 


In the above three methods, normally the maximum allowable limit for 
skewness is taken around 0.5 
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Aspect ratio 

This is the maximum ratio between of any two adjacent sides of a cell. Aspect 
ratio of less than 6 is desirable. It should not exceed 20. 

Adjacent cell ratio 

This IS a measure of grid smoothness. It is a ratio of the areas of two adjacent 
cells. It is desirable to have it close to one. 

Nearness Parameter 

This is a measure of how far the first grid lines are from the boundaries or how 
well does the grid intervals chosen at the boundary penetrate the interior. 

Fig 2.10 shows calculation of nearness parameter for a cell on the 7] = 
Vmax boundary. Let the distance between this boundary and its adjacent grid 
line be Sz at the cell and Sb at the corner. Also let the length of the grid line 
AB be Wi at the cell, and the length of the grid line at the edge CD be Wb- 
Then the nearness parameter is calculated as 

s^/w, 

Sb/Wb 


Nearness Parameter — 



2 6 Evaluation of Grid Quality 



Figure 2.10: Nearness parameter 


Among the two ^ = constant boundaries, ^ = ^rnm and ^ = ^max, the one 
which has the maximum St, is chosen for this calculation 



Chapter 3 


ELLIPTIC GRID 
GENERATION SYSTEMS 


In mid 1970s Thompson et al [15] first used elliptic Partial Differential Equa- 
tions for numerical generation of body-fitted coordinates. A number of grid 
generation packages were developed using Laplace or Poisson’s equations. Dif- 
ferent functions for controling grid point distribution in the interior were tried. 
These functions are essentially the source terms of the Poisson’s equations. 
Much effort was made by various authors to calculate the control functions 
automatically for achieving desired properties of the grid. 

In this chapter, elhptic grid generation systems and various control func- 
tions developed by various authors are discussed. 


3.1 Grid Generation Equations 

Elliptic grid generation is based on Poisson’s equations of the form 

(3.1) 

(3.2) 
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or Laplace equations of the form 


= 0 


V^77 = 0 


where 


= 


52 ^2 


+ 


dx^ dy 


(3.3) 

(3.4) 


By solving the above equations, using either Dirichlet conditions (grid points 
specified) or Neumann conditions (spacing specified) at the boundaries, we 
obtain as solution of the mapping 


^ = ^(x, y) , 


(3.5) 


and 


rj = r]{x,y) . (3.6) 

Once this is obtained we can find the inverse mapping 

X = x(^, 77 ) , (3.7) 

y = y{C,v) ■ (3.8) 

and find the (xp, yp) correspondence to each grid point P in the square trans- 

formed domain where (^p, rjp) are integer multiples of Af and A77, respectively. 
(A^ = 1/M ,A7] = l/A'") where M, N" are the number of grid intervals in the 
if and 77 directions respectively. Once we know (xp, yp) for each grid points, 
we have generated the grid and solved our problem. However , an easier way to 
do this is to convert Equations 3. 1-3.4 so that x, y are the dependent variables 
and f , 77 are the independent variables, as the values of the latter at the grid 
points are already known. On solving this transformed equations, we obtain 
Equations 3.7 and 3.8 directly, and thus know the grid in the physical domain. 
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So these equations are then transformed to the computational domain 
using the transformation relations discussed in the section 2.3, and the Equa- 
tions 3.1 and 3.2 transform to 

A{ax^^ - 2l3x^ri + + B{ay^^ - 2/?y^„ -f 7y„,) = P (3.9) 

C{<xx^i - 2/3x^r, + IXnn) + D{ay^^ - 20y^^ -f- 7y„,j) = Q 

(3.10) 

where 


A- 

B - 

ur 


C = 1 , 

D = 

1 J 3’ 

\J0 


0 = x^xr, -h y„y| 

7 = x| + y| 


and 

U 1= - ^vVi) 

Further simplification yields a system of two elliptic equations of the form 

- 20X^r, + JXr,^ = - \ J 0 (Px^ + Qxr) (3.11) 

- 20y^r, + IVvv = “ U T [PVi + QVr) (3.12) 

These equations are solved by finite difference methods on the uniform 
grid in (^, 77) space. These equations are non-linear since the coefficients 
a, /?, 7 are also functions of unknown variables x{^, 7), y(if, 7) and hence these 
equations need to be solved iteratively. Gauss Seidel iteration with under 
relaxation or successive line under relaxation can be used. We have used 
Gauss Seidel under relaxation in our program. 

In our computer program, all the second derivatives are centrally dif- 
ferenced. The coefficients a,/?, 7 are also calculated using central-difference 
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formulas using previous iteration values. To satisfy diagonal dominance of the 
discretised scheme, the first derivatives are differenced using upwind scheme 
i.e. if P is +ve then x^, are discretised using forward difference, else back- 
ward difference formula is used, and similarly x,,, y,, are discretised depending 
upon the values Q If the first derivatives are also approximated using central 
difference, the scheme will become unstable. 

The boundary conditions are given by initially choosing the grid points 
at the boundaries (Dirichlet conditions). Grid spacing and grid line slope can 
also be specified as the (Neumann) boundary conditions. But for internal flow 
applications, since all boundaries are fixed, only Dirichlet boundary conditions 
can be used. 


3.2 Properties of Elliptic Grid Generation Sys- 
tems 

3.2.1 Laplace Systems 

An important property of a Laplacian generated grid is its inherent smooth- 
ness. Due to this excellent smoothing property, boundary slope discontinuities, 
if any, are not propagated into the field. The Laplacian system also guarantees 
that the mapping obtained as the solution to the Elliptic equations is one-to- 
one. This property is very important to avoid ambiguity in the placement of 
grid points. 

With the Laplace system, grid lines will tend to be equally spaced in the 
absence of boundary curvature but it will become more closely spaced over 
convex boundaries and less closely spaced over concave boundaries. This lack 
of control over grid spacing may be problematic - often finer grid intervals are 
required in certain regions, eg. near solid boundaries etc., which laplace system 
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cannot deliver. This is the primary disadvantage that Laplace sysytems have 
in comparison with Poisson systems. 

3.2.2 Poisson Systems 

The source terms P and Q in the Poisson’s system can control the spacing 
between grid lines and the curvature of the grid lines at the boundaries. If we 
use Poisson’s equations with -ve values of Q near the boundary for the figure 
3.1, it will make rjyy more negative and hence r] — constant lines are further 
closely spaced near the boundary. Positive values of Q will have the opposite 
effort. 

Similarly the source term P controls the curvature (and hence the spac- 
ing) between ^ — constant lines. Since the boundary points are fixed, the 
— constant lines can not change its intersection points in the boundary. The 
overall effort of P is to rotate the ^ -constant lines in the direction of increasing 
^ at rj — constant boundary. 
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3.3 Control Functions 

Many methods are available to determine the control functions P, Q, in Pois- 
son systems. In this section some of them, which are used in our code, are 
discussed. 


3.3.1 Attraction of Grid Lines to Grid Lines/Points 


In the TOMCAT code of Thompson et al [12] they used source terms to attract 
the grid lines towards other grid lines/points. 


To achieve small spacing between ^ — constant lines in the neighbourhood 
of ^ = 1 ^ 1 , they used the control functions of the form 

n) = sign{i - ^,) exp {-B U - 6 1} 

The constant A is equal to the maximum magnitude of P and it controls the 
extent by which the ^ — constant lines shift towards if = ft while the range of 
the attraction effect is determined by the decay factor B. If a -ve value of A 
is used, instead of attraction, the grid lines will be repelled. 


To attract the ^-constant lines towards a given point (f., r},), the control 
function P takes the form 


^(f) 7 ) = - ft) exp |-5 [(f - ft)^ + (77 - ri,Y] ’ I 


(3.14) 


Similarly the functions 

77 ) = -A szgn{ri - r|^) exp {-B \ g - rjt \} (3-15) 

is used to attract g - constant lines towards 77 = 77 . line and the function 
Q(f, rj) = -A sign{g - 77,) exp [(f - f,)^ A (77 - 77,)^] ' | 


(3.16) 
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is used to attract 77 — constant lines towards the given point (if,, r;,). 

It should be noted that the constant A should not be too large, If it is, 
the mapping may no longer be one-to-one and grid lines of the same family 
(eg constant f lines) may cross each other. 

3.3.2 Thomas - Middlecoff Method 

Due to the strong smoothing effect of the Laplacian operator the grid lines in 
the interior will tend to be equally spaced irrespective of the boundary point 
distribution. This is not always desirable. In 1980, Thomas -and Middlecoff 
[11] proposed a direct control method whereby the resulting interior grid point 
distribution is controlled entirely by a prior selection of the grid point distribu- 
tion along all the boundaries. The angle of intersection between the transverse 
grid lines and the boundaries can also be controlled. 

Thomas and Middlecoff chose source terms of the form 

p = <!>((. r,) (d + e,) 

Q = vCf. ij) i’ll + nl) 

where the functions (p, ip will be calculated from boundary point distribution. 
Introducing these terms, the Poisson’s Equations 3.11 and 3.12 assume the 
form 


q:(x 4^ + (px^) - 2/3x^^ + 7(x,,„ + ipXr,) = 0 (3.17) 

Q:(y^{ + + liVnv + ^Vn) = 0 ( 3 - 18 ) 


The functions (p and ip are obtained by imposing two constraints i.e. the grid 
lines transverse to the boundary axe (a) locally straight and (b) orthogonal to 
the boundary. 
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To find an expression for (f) along rj = rjb = constant boundaries, is 
eliminated from the Equations 3.17 and 3.18. This yields 


[Vv 


+ <py€)] = vl 




(3.19) 


The ratio ^ represents the slope of the ^ — constant lines that are transverse 
to the boundary r] = 775 . The first constraints is that the transverse coordi- 
nate lines is straight (i.e. have zero curvature) in the neighbourhood of the 
boundary. Mathematically 


dr) 




= 0 at T) = r)b 


The second condition is that the transverse grid lines are locally orthogonal to 
the boundary 77 = 775 . Let a vector r = (x, y) denote a point (x, y) in physical 
domain. Then the vector that is locally tangent to 77 — constant line is 


Ti = Vi) . 

Similarly the local tangent vector to ^ — constant line is 

— (^ 17 ? J/r?) • 

The two families of coordinate lines are orthogonal if and only if 


7-^.r„ = 0 


or 

+ VcHv = 0 (3.20) 

Under this orthogonality condition 0 = 0. Hence both the terms in the right 
hand side of the Equation 3.19 vanish at the boundary 77 = 774 . The Equation 
3.20 is used to eliminate the first derivatives x,,, y,, from the Equation 3.19. 
This yields an expression for 0 as 

, - + y^y«) 


at r) — Vb 


(3.21) 
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Similarly i/; is obtained as 

[^n + yn) 

The values of 4>, t/’ in the interior are obtained by linearly interpolating 
</} in rj direction and ^ in ^ direction. All the first and second derivatives in 
3.21 and 3.22 are approximated using central difference. 

The evaluation of the functions ip is equivalent to constructing a local 
curve-fit to the arc length between the pre-assigned boundary points. The 
interpolation extends the range of curve-fit. Anderson and Steinbrenner [1] 
proved that the Thomas-Middlecoff formulation does not completely reflect 
the curvature effect of the boundary for interior points near boundaries. How- 
ever, the Thomas-Middlecoff method gives satisfactory results for most of the 
problems. Also it converges quickly compared to Sorenson and Steger method 
which will be discussed in the next section. If the boundary points are placed 
at equal intervals, then the grid generated by Thomas-Middlecoff method and 
Laplace equations are same. 

3.3.3 Sorenson Method 

For some applications it is desirable to have grids that are orthogonal at the 
boundaries, for then the boundary conditions reduce to simple equations and 
certain terms may vanish. With nonorthogonal grids, the boundary conditions 
may transform into cumbersome equations in the computational domain, and 
may induce loss of numerical precision. 

In 1979, Steger and Sorenson [8] developed a method to calculate the 
source terms automatically so that the intersection angle between the trans- 
verse coordinate lines with the boundaries, and the spacing between the bound- 
aries and their adjacent coordinate, line can be controled effectively . This 
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essentially needs a judicious choice of the source functions P and Q in Equa- 
tions 3 11 and 3.12. In this method, the spacing between the boundary and 
the adjacent grid line is specified. Also specified is the angle between the 
boundary line and its transverse lines. Using these constraints, the values of P 
and Q at the boundaries are determined by a prescribed procedure. Knowing 
these values at the boundaries, the P, Q values in the interior domain can be 
determined by interpolation. After that, the Equations 3 11 and 3.12 are used 
to generate the grid. 



Figure 3.3: Sorenson Method 


Suppose, the spacing and the intersection angle are to be controlled at 
rj — constant boundaries (77 = rimm and rj = rimax)- The P and Q values in the 
interior are chosen as 


P T}) = Pr,„..„ (0 exp {-a (r? - Vrmn)} + Prima. (0 exp {-b iVmax -V)} 

Q (^, v) = Qvm,n iO exp {-c (77 - Vmzn)} + Qvma. iO exp {-d {rimax - v)} 

(3.24) 


where a, b, c, d are some positive constants, whose value has to be spec- 


ified. 
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The first constraint imposed on the grid system is that the spacing along 
^ — constant lines between the boundary rj =: rji, and its adjacent grid line is 
specified. Let this desired spacing be As then 

As = l^(Ai)^ + (Aj/)^j ^ at rj = T}i, 

where Ax and Ay are the differences in x and y over the interval. When Ax 
and Ay approach zero, the differential relationship is 

ds = [(dx)^ + {dyf\ ^ at rj = 

Using the transformation relations 

ds = [{x^d^ + Xr,dyf + (y^d^ + yr,dT]f ] " at V = Vb 

Since is constant along the interval under consideration, the above equation 
reduces to 

ds = (xl + y^) ^ drj at J? = 77 ;, . 

or equivalently 

5 .= + = ( 3 - 25 ) 

Sr, is specified along the boundary r] = rjb. 

The second constraint imposed on the grid system is that the angle of 
intersection between the transverse grid lines and the boundary is specified. 
Let the angle be 9 \n=nb- Then from dot product definition 

V^V 77 = I Vf I I Vy I cose 

Expanding this equation and transforming it to the computational domain we 
get 

-^vH-ynVi = + {xl + ylY cos 9 


(3.26) 
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Similarly from the cross product definition 


Vif X V 77 = I Vf I I 7 I smd 


+ yj) ^ (a:| + yf) " sin 6 


(3.27) 


Solving Equations 3.26 and 3.27 for and j/^ and using Equation 3.25 we get 

^ ^ ^ 23 ) 

N + !/|)" 


(3.28) 


(— y^cos^ + x^sinO) 


^l + yp^ 


at Tj = r]b . 


(3.29) 


Both 9 and s, is specified at r; = 77 ;, and can be calculated by 

differencing the x, y coordinates of the chosen boundary points. So x^, y„ can 
be calculated zX vj = rji from the above two equations. 

With similar procedure on ^ = ih boundaries, we can get x^, y^ as 


x^ 


■Xr, COS 6 + y,,sin^) 

at 


II 

(3.30) 



(^n + VvY 







■yr, cos 6 — x,,sin^) 

at 


= 6 - 

(3,31) 


+ y^) 


where + y|) ^ is the spacing between ^ boundary and its 

adjacent grid line which needs to be specified. 

The four positive constants a, b, c, d in the Equation 3.23 and 3.24 
determine the extend by which the control functions at the boundaries affect 
the interior. These constants are chosen such that the effect of control functions 
of the one boundary is vanishingly small at the opposite boundary, so only one 
term will remain in the right hand side of the Equations 3.23 and 3.24 along 



3.3 Control Functions 


34 


Tj — r]b boundaries. Now substituting into the Equations 3.11 and 3.12 we get 


where 


V) 

Q{^, v) 


VtiRi ^17-^2 

[T\ 

\T\ 


= 

R2 = 


\j? 

- (ay« - + 7y.w) 

UP 


at rj — r}b (3.32) 

at t? = % (3.33) 

at T] = T]b (3.34) 

at r} = rji (3.35) 


We will obtain the same Equations 3.33 - 3.35 at ^ if the grid lines 
were to be controlled at ^ 


For boundaries rj = rjb, The derivatives x^, y^, are simply 

calculated from central difference approximation of the chosen boundary points 
and the first derivatives x^ and axe calculated using the Equations 3.28 and 
3.29. The cross derivatives x^„, y^ri can be computed by differencing x„ and y,,. 
The double derivatives x,,^ and y,,,, are computed by differencing the current 
approximation of the x, y field using 


_ -7n,., + 8x|,., - i|„, _ 

“ 2 (A,)" A, 

+ 8y|,., - yi,., _ 

“ 2(A,)= At, 


~'^^\j= 3 max ^|j=jmoJ-2 ^ ^|j=j max 


+ 


2{Ayf 

_ ~f 8y|j-j^oj_i — y — |j=jmai -2 


y‘^'^\j=3max 2 {Arj)^ 

where the grid node z, j corresponds to 


Ay 


+ 


Ay 


!(•'. j) = (•-!) Ae 

7 (■. j) = y - 1) 
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Because x^r, and need t he knowledge of the x, y field, this method needs an 
iterative solution procedure. The double derivatives a:,,,, and yj^ri keep changing 
with iterations of the the solution until convergence is achieved. 

The source terms can be quite large and if in every iteration, values 
calculated from Equations 3.33 and 3.33 are used exactly, instabilities can 
result. So under relaxation becomes essential in the iterative procedure. The 
source terms at t he boundary are updated by using the following expressions. 

pin+l) ^ pin) ^ Jj-1 _ />(«)' 

g(n+l) ^ Qin) ^ ^ - gW] 

where n is the iteration level and the weight functions wp, wq are quite small 
(0.02 to 0.06). 

For sharp corners, occuring along the boundaries, the computed source 
terms have to be replaced by an average of the source terms on either side of 
the corner. 

The above method, as discussed by Steger and Sorenson [8], can control 
the grid lines over only two boundaries. To control the grid lines over all the 
boundaries, instead of control functions of the form Equation 3.23 and 3.24, 
we used 

r]) = Pir^.Av) T? (1-0(1-^) exp{-a(^-^,„.„)} 

+ Pua. iv) Cl (1 - l) exp {-6 {Uax - 
+ Prirmn iC) ^ (1 “ 0 (1 “ l) G^V{-<^{l “ 1mA} (3.36) 
+ PvmaxiC) Cl{^~C) exp {— d (77 toox ~ ^)} 

Q iC, l) = Qu,n il) ^ (1 - 0 (1 - l) exp {-a iC - Un)} 

+ Quax il) Cl{^-l) exp {-b {(max - 0} 

+ Qvmm iO e (1 - 0 (1 - l) exp {-c{l- lm,n)} (3.37) 

+ QvmaxiC) Cl {1-0 exp{-d{7)max-l)} 
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In many geometries, it is impossible to maintain both orthogonality and 
small spacing between grid lines in all the boundaries. If we use the Sorenson 
method for these geometries, the solution will not converge or a bad grid will 
result even after a thousands of iterations. 

In our experience, we found that instead of calculating both source terms 
P and Q in all the boundaries, calculating P in ^ — constant boundaries and 
Q m Tj — constant boundaries yields good results. Let us call this method as 
Modified Sorenson method. The evaluation of source terms at the boundaries 
is the same as that of the Sorenson method, except that only one source term 
is determined from each boundary. The control functions in the interior are 
calculated from 

PU,V) = Pirmniv) exp{-a(^-^m.n)} + Pua.iv) 

= Qrtm.AO exp{-c(r/-r?,„.„)} + QnmasiO 


As we know that the source term P controls the spacing in the ^ direction 
and intersection angles at the rj- constant boundaries and Q controls the spac- 
ing in the ry direction and intersection angles at the ^ - constant boundaries, 
the Modified Sorenson method stresses more on spacing control than control of 
intersection angles. Even so the Modified Sorenson method generally delivers 
a near-orthogonal grid in the boundaries. 

In the Modified Sorenson method, to save computations, the control 
functions are calculated only once instead of recomputing x,, in r, - constant 
boundaries and in ^ - constant boundaries are updated after each Gauss 

seidel iteration sweep 
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3.4 Orthogonal Grids - Moving Boundary Points 
Method 


Without using any control functions, orthogonality can be achieved at the 
boundaries by suitably choosing the boundary points. But there is no direct 
way of choosing the correct grid points without knowing the interior grid points 
which themselves depend upon the boundary grid point distribution. So the 
boundary points are chosen, and then corrected iteratively to obtain a locally 
orthogonal grid at the boundaries. This method, given below, is called the 
Moving Boundary Points Method. 

To be orthogonal, the grid should satisfy 

x^Xn + ViVT) = 0 (3.40) 

at that point. First with some assumed boundary point distribution, interior 
grid points are calculated by solving the Laplace Equations 3.3 and 3.4. To 
find out new boundary points on ^ — constant boundaries, if | x,, [ < | j/,, I 
then is calculated from the orthogonality condition 3.40 i.e. 

where the quantities on the right are estimated from the previous assumed 
X, y field. If I yj < I xj, then x^ is similarly computed 

ViVv 

x^ • 

Xjj 

This procedure avoids the difficulty caused if either x„ or yr, is zero. Suppose, if 
is calculated, then the x coordinate of the new boundary point is recomputed 
from the finite difference form of x^ keeping the interior grid nodes at previous 
values. Then y is computed from the equations of the boundary. Suppose if 
% is calculated first, then y is calculated from the difference from of and x 
is calculated from the equations of the boundaries. 
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It is necessary to slow down the corrections in the boundary points near 
the corners. Otherwise some of the points in one boundary can move to the 
other boundary and the mapping will not be one to one anymore. So an under 
relaxation factor which is a function of the distance of the point from its nearest 
corner is used. Thus the new coordinate is calculated using 


^new 




where 

pre ■ represents previous value. 

cal - represents the value calculated from the finite difference form. 
d - is the distance between the grid point and its nearest corner. 
dn - can be chosen as one third of the distance between two nearest 
k - is some constant in the range 1-4. 


The y coordinate corresponding to the calculated x is calculated from the 
equations of the boundary. 


For rj — constant boundaries, a similar procedure is used. This whole 
procedure is repeated until the desired grid is obtained. This method is compu- 
tationally very expensive compared to other methods. We have implemented 
this method only for a particular section. 



Chapter 4 


RESULTS AND DISCUSSION 


While writing codes for grid generation, we are most concerned about bound- 
ary orthogonal grids which will make implementation of boundary conditions 
simple and in the clustering of grid points near the boundaries to capture the 
boundary layer. However there are many other considerations which need to 
be included to certify a grid a.s ‘good’. It is difficult to identify a single method 
which is the best of all because the requirements of the grids vary from problem 
to problem. 

Grids generated by various methods for the same domain can be com- 
pared using the parameters such as skewness, aspect ratio, adjacent cell ratio, 
nearness parameter, etc. Even then if the parameters are within certain lim- 
its, all the grids are acceptable: at the same time there is no unique general 
criteria with which to certify a particular grid the best. 

A method which gives a better value for a particular property may not 
give a better result for a second property . In such cases, a greater relative 
importance can be given to certain parameters depending upon the problem 
in which the grid is to be employed. For examples, for viscous flow problems, 
the nearness parameter can be very important since it ensures relatively more 
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number of grid points near the boundaries. 

In all the methods, the interior grid point distribution depends upon 
boundary grid point distribution which is entirely specified by the user. So the 
intitutive skill in the choice of boundary points also comes in to play. 

The quality of the grid generated by the Sorenson and Modified Sorenson 
method depends upon the exponential decay factor, the spacing between the 
boundaries and their adjacent grid lines, and the intersection angle between 
the grid lines at the boundaries, which are also specified by the user. 

To study the effectiveness of various methods, grids were developed for 
the following sections. 

1. Circular section 

2. Converging section 

3. L - section 

4. Quarter portion of a super ellipse 

5. Channel section with a downward slope at the middle 

Grids were developed for both uniform as well as nonuniform boundary point 
distribution. Except for the circular section, exponential functions were chosen 
to generate nonuniform boundary point distribution. The various parameters 
viz skewness, aspect ratio, adjacent cell ratio, nearness parameter are com- 
puted for all the grids and were used for comparison. These parameters are 
defined in section 2.6 of this thesis. 
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4.1 Case 1: Circular Section 

The figures 4.1(a)-(e) show grids for a circular section with equally spaced 
boundary points. The figures show the grids obtained by (a) Transfinite in- 
terpolation (b) Laplace method (c) Thomas-Middlecoff method (d) Sorenson 
and (e) Modified Sorenson method, respectively. The salient features of the 
various methods emerge from a comparison of these figures. Transfinite inter- 
polation gives highly skewed cells at the corners, The Laplace grid (Fig 4.1(b)) 
has equally spaced and orthogonal grid lines except near the boundaries from 
where the grid lines were repelled due to concave curvature of the boundary. 
With equally spaced boundary points, Thomas-Middlecoff method 4.1(c) gives 
the same results as that of the Laplace system. The Sorenson method 4.1(d) 
gives a near-orthogonal grids in the boundaries. Modified Sorenson method 
tries to get the specified spacing at the boundaries. For both Sorenson and 
Modified Sorenson method, half of the spacing between the boundaries and 
their adjacent grid lines of the Transfinite grid is specified as the spacing. Also 
the intersection angle at the boundary is chosen as 90°. The exponential decay 
factor is taken as 0.5/ Arj. 

Figures 4.2 (a)-(e) shows the plot of cumulative percentage of the grid 
cells Vs various parameters. All the parameters - skewness, adjacent cell ratio, 
aspect ratio, intersection angles between grid lines, nearness parameter - were 
taken at the abscissa. The plots show what percentage of the grid cells are 
above or below that value on the horizontal axis. Adjacent cell ratio and aspect 
ratio are nearly same for all methods. The Laplace method and Thomas- 
Middlecoff give less skewed grids while Sorenson gives the most skewed grid. 
This is due to the cells near the corners. The Modified Sorenson performance 
regarding skewness is somewhat better than that of the Sorenson method. The 
nearness parameter is nearly same for all methods. 
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Figure 4.3(b) shows the grids generated by iteratively moving the bound- 
ary points to get boundary orthogonal grids, as discussed in section 3.4. The 
boundary points were moved towards the corners. The figures 4.3(a) and 4.3(c) 
show the grids generated by the Transfinite interpolation and the Thomas- 
Middlecoff method for the same boundary points as in figure 4.3(b). Both 
Sorenson and Modified Sorenson methods give bad grids for this boundary 
point distribution. For these results, the parameters were plotted in the fig- 
ures 4.4(a)-(e). 

From these figures, it can be seen that Transfinite interpolation gives 
highly skewed grids at the corners. Adjacent cell ratio is minimum for the 
Laplace method which is due to its strong smoothing property. However, 
although the Laplace method gives less skewed, smooth grids, its grid lines are 
far away from the boundary. Hence Laplace method is not suitable for viscous 
flow calculations. The Thomas-Middlecoff gives a fine grid at the boundaries, 
suitable for resolving the boundary layer. This is undoubtably the best grid. 


4.2 Case 2: Converging Section 

The figures 4 . 5 (a) and 4 . 7 (e) show grids generated for a converging section 
with two types of boundary point distribution. The figures 4.6 and 4.8 show 
the properties plots for this grids. With equally spaced boundary points, all 
methods are doing well. For both Sorenson and Modified Sorenson method, 
the exponential decay factor is chosen as I/A 77 . and the spacing is specified as 
0.75 times the spacing in the Transfinite grid. It can be easily seen from the 
figures that due to higher value of exponential decay factor only few grid lines 
near the boundary were affected. Nearness parameter is equal to 1 for the grid 
generated by Transfinite interpolation. For other methods, it varies from 0.5 

to 1.5. 
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ThomcLS-MiddlecofF is doing well with exponentially spaced boundary 
Thomas-Middlecoff gives the results with closely spaced near-orthogonal grid 
lines at the boundaries. Transfinite interpolation itself gives a better grid than 
the Laplace, Sorenson and Modified Sorenson, which all tend to have large 
spacing near the boundaries. 


4.3 Case 3: L - Section 

The figures 4.9(a)-(e) and 4.11(a)-(e) show the grids generated for a right 
angled bend. The boundary discontinuity was propagated into the grids gen- 
erated by Transfinite interpolation. For the equally spaced boundary points 
case, the Transfinite interpolation seems to be best for the others gives rather 
large grid spacing near the corners. But on the whole, the grids are compared 
in terms of skewness, aspect ratio etc.. In the exponentially spaced boundary 
points case, Thomas-Middlecoff does best in terms of the grid performance 
near the boundary with Transfinite interpolation coming second, the others 
showing unacceptable performance at the corners. 

The figures 4.10 and 4.12 show the property plot for these grids. Trans- 
finite interpolation gives grids with lower skewness. Sorenson method gives 
higher number of grid cells with slightly large skewness. Adjacent cell ratio 
is slightly better for Transfinite grids. With equally spaced boundary points, 
Transfinite and Laplace method gives smaller aspect ratio grids than Sorenson 
or Modified Sorenson method. But with exponential boundary point distri- 
bution, Thomas-Middlecoff and Transfinite interpolation give large number of 
grid ceUs with higher aspect ratio, but from figures 4.10(a)-(e) they give the 
best grids!. All this goes to show that it is difficult to judge the performance 
of the methods by quantitative parameters - the best judge of god grid may 
still be qualitative ‘eye ball’ estimate. 
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4.4 Case 4: Rectangular Channel Section with 
a Slope 

The figures 4.13(a)-(e) and 4.15(a)-(e) show the grid generated for a channel 
with a downward slope at the middle. As before, we can see the discontinuity 
in the interior in the grid generated by Transfinite interpolation. Also the grid 
lines are away from the boundary near the concave portion of the boundary in 
the grid generated by Laplace method. Sorenson method tries to get orthogo- 
nal grid lines and hence make the grid lines curved which are vertical in other 
methods. 

Figures 4.15(a)-(e) and 4.16(a)-(e) show the case of exponentially spaced 
boundary points. As usual, with exponentially spaced boundary points, Thomas- 
Middlecoff method does well. Transfinite interpolation also yields an accept- 
able grid except in terms of smoothness. The other methods fail in maintaining 
nearness at the boundary. However both Thomas-Middlecoff and Transfinite 
interpolation give rather large aspect ratios (figure 4.16(c)). 

4.5 Case 5: Quarter Portion of a Super Ellipse 

The figures 4.17(a)-(e) and 4.19(a)-(e) show the grids generated in a quarter 
portion of the superellipse (f) ~ Iwithn = 3, ^ = 1.4. The same 

illustration was used by Thomas and Middlecoff [11] to test their method. 

With equally spaced grid points in the boundaries, Transfinite, Laplace, 
Thomas-Middlecoff and Sorenson give similar results. While Modified Soren- 
son fails to give a good grid. With exponentially spaced boundary points 
Thomas-Middlecoff gives the best grid with evenly spaced cells close to the 
boundary. Transfinite interpolation also gives an acceptable grid while the 
Laplace tends to pull the grid lines away from the boundary. 
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4.6 Conclusions 

A number of methods of grid generation, Transfinite interpolation, Laplace, 
Thomas-Middlecoff, Sorenson, Modified Sorenson method, were used to gen- 
erate H-type grid for a variety of configurations and their performance was 
evaluated. Special attention was paid to the ability of the metiiod to deliver 
grid lines that followed the boundaries closely. This, property has great rele- 
vance for grid generated for viscous computations in internal geometries. The 
performance of Thomas-Middlecoff was clearly the best with Transfinite in- 
terpolation coming second, although the grids of the later were not always 
smooth. The other methods were quite erratic, often failing badly. 
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Figure 4.1; Circular Section with Uniform Boundary Point Distribution 
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Figure 4.2: Parameters Plots for circular Section with Uniform Boundary Point 
Distribution 
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(c) Thomas-Middlecoff Method (d) Sorenson Method 



(e) Modified Sorenson Method 
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4.3: Circular Section with Nonuniform Boundary Point Distribution 
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Figure 4.4: Parameters Plot for Circular Section with Nonuniform Boundary 
Point Distribution 








Cummulalive % of grid cell Cummulative % of grid cell 


4 6 Conclusions 






Nearness Parameter 


(e) 

Figure 4.6; Parameters Plot for Converging Section with Uniform Boundary 
Point Distribution 









Cummulatjve % of grid cel! Cummulalive % of gnd cell 


4.6 Conclusions 


53 




(a) 


(b) 




(c) 


(d) 



(e) 


Figure 4.8; Parameters Plot for Converging Section with Nonuniform Bound- 
ary Point Distribution 







4 6 Conclusions 


54 



(a) Transfinite interpolation 



(b) Laplace method 




Figure 4.9; L - Section with Uniform Boundary Point Distribution 
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Figure 4.10: Parameters Plot for L - Section with Uniform Boundary Point 
Distribution 
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(c) Thomas-Middlecoff Method (d) Sorenson Method 



(e) Modified Sorenson Method 


Figure 4.11: L - Section with Nonuniform Boundary Point Distribution 
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Figure 4.12; Parameters 
Distribution 
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(c) Thomas-Middlecoff Method (d) Sorenson Method 



(e) Modified Sorenson Method 


Figure 4.13: Rectangular Channel Section with a Slope with Uniform Bound- 
ary Point Distribution 
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Figure 4.14: Parameters Plot for Rectangular Channel with a Slope with Uni- 
form Boundary Point Distribution 







(a) Transfinite interpolation 


(b) Laplace method 



(c) Thomas-Middlecoff Method 


(d) Sorenson Method 



(e) Modified Sorenson Method 


Figure 4.15: Rectangular Channel with a Slope with Nonuniform Boundary 
Point Distribution 
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Figure 4.16: Parameters Plot for Rectangular Channel Section with a Slope 
with Nonuniform Boundary Point Distribution 
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Figure 4.18: Parameters Plot for Quarter Portion of a Super Ellipse with 
Uniform Boundary Point Distribution 
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Figure 4.19: Quarter Portion of a Super Ellipse with Nonuniform Boundary 
Point Distribution 
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(c) 
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Figure 4.20: Parameters Plot for Quarter Portion of a Super Ellipse with 
Nonuniform Boundary Point Distribution 
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